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Abstract. We present a new type of search strategy designed specifically to find 
continuously emitting gravitational wave sources in known binary systems based on the 
incoherent sum of frequency modulated binary signal sidebands. The search pipeline 
can be divided into three stages: the first is a wide bandwidth, JF-statistic search 
demodulated for sky position. This is followed by a fast second stage in which areas 
in frequency space are identified as signal candidates through the frequency domain 
convolution of the ^-"-statistic with an approximate signal template. For this second 
stage only precise information on the orbit period and approximate information on 
the orbital semi-major axis are required apriori. For the final stage we propose a fully 
coherent Markov chain monte carlo based follow up search on the frequency subspace 
defined by the candidates identified by the second stage. This search is particularly 
suited to the low-mass X-ray binaries, for which orbital period and sky position are 
typically well known and additional orbital parameters and neutron star spin frequency 
are not. We note that for the accreting X-ray millisecond pulsars, for which spin 
frequency and orbital parameters are well known, the second stage can be omitted and 
the fully coherent search stage can be performed. We describe the search pipeline with 
respect to its application to a simplified phase model and derive the corresponding 
sensitivity of the search. 



PACS numbers: 04.80.Nm, 07.05.Kf, 95.55.Ym, 97.60.Gb,97.80.-d 



1. Introduction 



The LIGO and GEO600 ground based interferometric gravitational wave detectors have 
been taking data in "Science runs" of ever-increasing length and sensitivity, since August 
2002. At present the detectors are midway through S5, a year long run started in 
November 2006, with LIGO running at design sensitivity and GEO600 not far behind. 
Despite the unprecedented sensitivity of these instruments the data analysis task of 
identifying weak continuous gravitational wave (GW) signals buried in the detector 
noise represents a vast computational challenge. To date a number of analyses have been 
performed using both LIGO and GEO600 data searching for targeted sources such as the 
known isolated and binary radio pulsars [U [2] , and for wide parameter space searches 
including all-sky searches for isolated rapidly rotating neutron stars [3l H] and searches 
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for GWs from the low-mass X-ray binary (LMXB) Sco X-l [U [5]. Present matched 
filter based search techniques, although optimal in terms of sensitivity for a given 
integration time, are unfeasibly computationally intensive due to the steep rise in the 
number of search templates with observation time. As such, new less computationally 
intensive sub-optimal incoherent and hierarchical techniques must be developed whereby 
sacrifices in signal-to-noise ratio (SNR) can be recovered through longer observation 
times. Examples of such methods include the stack-slide [6] , Hough [3] , and radiometer 
[5] searches. 

We present a method previously employed in the electromagnetic detection of radio 
pulsars [7J which we have adapted specifically for the detection of continuous GW sources 
in binary systems such as LMXBs. This method, which we call the "sideband" search, 
is based on the incoherent summation of signal power present in the finite number 
of frequency modulated (FM) sidebands expected from such sources. This "sideband" 
search stage is used to identify signal candidate frequencies forming a reduced frequency 
subspace on which a coherent Markov chain monte carlo (MCMC) based follow up can 
be run. We identify that this follow up stage alone is suitable for searching for the 
accreting millisecond X-ray pulsars (AMXPs) [SI El EH QH H21 031 03! for which the 
frequency is already known. 

We should stress that for the purposes of these proceedings we provide a description 
of this method applicable to the very basic case of circular orbit monochromatic sources 
of exactly known orbital period and sky position and assume a constant and flat detector 
noise spectrum within the band. Here we aim to establish the principles behind this 
search strategy but for more general cases details of the practical application of this 
pipeline can by found in [T5], [16] . 



2. GW emission from LMXBs 



Observations by the Rossi X-ray Timing Explorer (RXTE) have provided evidence 
supporting the idea the that accretion torque supplied to neutron stars in accreting 
binaries could provide a sustained source of angular momentum able to power GW 
emission. By balancing the accretion torque with the angular momentum radiated from 
the system through GWs we can estimate the GW amplitude received at Earth as 

" ~ 8 5 X fe) 1 (look) ' ( lO-Lr- ) 2 ' (1) 

where r is the distance to the binary, / is the gravitational wave frequency (twice the 
rotation frequency) and M is the accretion rate inferred from observed X-ray luminosity. 
The RXTE observations showed an apparent clustering of spin frequencies in the range 
250—370 Hz well below the theoretical neutron star break-up frequency of ~ 1 kHz. This 
implies the necessity for a mechanism through which an equilibrium can be achieved 
with the accretion torque. The competing (non GW based) explanation for this possible 
equilibrium state is provided by an accretion disk - magnetosphere interaction model, 
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[T7] in which there is required a relation between accretion rate and neutron star 
magnetic field strength. As this relation is not expected, the argument for GW emission 
from these systems remained plausible [H] and has motivated the detailed investigation 
of the ability of a neutron star to support the non-zero quadrupole moment (or current 
quadrapole in the case of "r- modes" [19]) needed for GW emission [20J. More recently 
work by [21] using a refined accretion model has been able to remove the need for an 
additional spin-down torque and although this therefore removes the need for GWs from 
these sources to explain the observed spin rates it does not make the GW generation 
mechanisms any less viable. 

LMXBs are a class of semi-detached binary system consisting of a neutron star or 
black hole in orbit around a lower mass Roche lobe filling companion object, usually a 
white dwarf or brown dwarf star. They have long been thought to be strong candidates 
for continuous gravitational wave emission due to the large (~ 100% Eddington limit in 
some cases) accretion rates inferred from X-ray luminosity. There are 85 know LMXBs, 
[22], with periods ranging from ~ 700 sec to ~ 33 days and of these, 10 have been 
seen to exhibit type 1 X-ray bursts [23] allowing us to measure the spin frequency of 
these sources to an accuracy of ~ 1 Hz. It has long been thought that the separation 
frequency between pairs of the kHz quasi-periodic oscillations (QPOs) seen in 17 of the 
known LMXBs are directly related to the neutron star spin frequency [24J and recent 
observations of kHz QPOs in the millisecond accreting X-Ray pulsars SAX J1808.4- 
3658 [25] and XTE J1807-294 [26] have lent weight to this theory. Unfortunately in 
LMXBs the kHz QPO separation frequencies are not constant and are seen to vary with 
source brightness. Indeed for the 7 LMXBs for which bursts and kHz QPO separation 
frequencies are known the simple relations obeyed by the AMXPs is not so clear cut. 
As a consequence we should assume that the unknown LMXB spin frequencies obey 
the same distribution described by those systems for which the frequency is known, 
currently ranging 250 — 620 Hz. This implies that any sensible and exhaustive search 
for these objects should use the entire ground-based interferometer sensitivity window 
of ~ 100 - 1500 Hz. 

3. The signal model 

We begin by assuming that the data received at a GW detector can be represented as 
the sum of the signal and additive Gaussian noise such that 

s(t) = W(t)[h(t)+n(t)\, (2) 

where the function W(t) is a time domain window function required to describe the 
non-continuous operation of a gravitational wave detector and can be equal only to zero 
or unity at any given time. We consider GWs emitted via quadrupole radiation from 
a rotating non-axisymmetric triaxial ellipsoid (GWs are emitted at twice the rotation 
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frequency). Following [27] we write the signal as 

4 

h(t) = Y,AiH*), (3) 

i=l 

where the signal amplitudes A\ are 

Al — [ho (l + cos 2 l) cos 2ip cos $o — 2/i cos t sin 2?/> sin $ ] /2, 

^2 = [/io (l + cos 2 i) sin 2-0 cos $o + 2/io cos i cos 2-0 sin $ ] /2, 

^3 = [ — (l + cos 2 t) cos 2^ sin $ — 2/i cos i sin 2^ cos $o] /2, 

At = [— /io (l + cos 2 i) sin 2-0 sin $ + 2/i cos 6 cos 2?/> cos $ ] /2. (4) 

Here /iq is t ne signal amplitude, t is the inclination angle (the angle between the NS 
spin axis and the line of sight vector from source to detector), if> is the GW polarisation 
angle, and <£> is the signal phase at t — t . The 4 corresponding time dependent signal 
components h(t)i are 

h\ — a(t) cos$(t), h 3 = a(t) sin , 

fr 2 = 6(t)cos$(t), /i 4 = &(t)mn$(t). (5) 

where $(t) is the time dependent GW phase and the functions a(t) and b(t) are related 
to the GW detector antenna pattern functions F + and F x by 

F + = a(t) cos'0 + b(t) simp, F x = b(t) cos ip — a(t) sin ip. (6) 

3.1. A simple binary phase model 

We have limited our phase model to that of a monochromatic source in a circular binary 
orbit, so the corresponding GW phase is 

= 27if {t - t ) + 27cf a cos (2vr(t - t p )/P) , (7) 

where fo is the intrinsic (constant) GW frequency of emission, a is the light crossing 
time of the orbital radius projected along the line of sight, P is the orbital period, and 
t p represents a reference time at which the source passes through the ascending node 
of the orbitQ- We consider known systems for which the sky position is known to high 
accuracy and as such we assume that any phase contribution from the detector's motion 
with respect to the solar system barycentre (SSB) can be removed. 

3.2. The Jacobi-Anger identity 

The chosen phase model in equation [7] shows that we are dealing with an FM signal, 
with modulation amplitude 27if a, and period equal to the orbital period P. Using the 
Jacobi-Anger identity 

oo 

e izcosd = J2 i n Jn{z> [nd i (8) 

n=— oo 

| For eccentric orbits t p represents the time of orbital periapsis, hence the subscript p. 
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we can express this as an approximate complex phase factor 

e^) » f; J n (Z) exp { 27 rit (/„ + £) + in (| - |) } , (9) 

n=—m *■ \ / J 

where J n is the Bessel function of the first kind with argument Z = 2nf Q a. Note that we 
have reduced the problem from an infinite to a finite summation by taking advantage 
of the fact that J n (Z) « for n > Z. In practise our approximation will be set to 
m = Z = luf^a. 



3.3. The T -statistic 



The jF-statistic [27] is obtained through maximisation of the log likelihood with respect 
to the signal amplitudes Ai. It is constructed from the complex components 

F a {f)= r a^h^W^e-^dt, (10) 



F b (f) = [ 3 bitM^W^e-^dt, (11) 
Jo 

where the total observation time span is T s and the total observation time T Q is 
Iq Ts W(t)dt. The JF-statisti(J§| itself is then defined as 

mf) = To sl U)D [ B \ F M)\ 2 + A\F b (f)\ 2 - (F a (f)F b (fT)} , (12) 
where the quantities A,B,C,D are 

A = |- r a 2 (t)W(t)dt, C = |- P a{t)b(t)W(t)dt, (13) 

1 o Jo 1 o Jo 



5 = J- / b'\t)W(t)dt, D = AB- C'\ (14) 



T 

- 1 - O 

and Sh(f) is the single-sided noise spectral density. By applying our phase model 
(equation |7J) we can write the components of the ^-"-statistic as 

F a (f) =(A + C) (A 2 - L4 4 ) jr ^e^-^Wif - f n ) , (15) 

n=—m 

F b (f) =(C + B) (A, - L4 3 ) f; ~*)w (f - f n ) , (16) 

n=—m 

where f n = /o + n/P and we have neglected negative frequency components. We have 
used the approximations 

P iy(t)a 2 (t)e" 2 ^dt « ]-AT W{f), (17) 
/ ° W^it^-^dt « isr o W(/), (18) 
y ° ^(t)a(t)6(t)e" 2 ^dt « ^CT W(/), (19) 



Here we choose to use IT rather than T which is statistically more easily described. 
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Figure 1. The jF-statistic as a function of frequency for a simulated signal with 
parameters f = 200 + 1.0849 x 10~ 6 Hz, a = 0.005 sec, P = 7912.7 sec, a = 2.0 rad, 
5 = 1.0 rad, and nuisance parameters ho = 1, if> = 0.2 rad, cost = 0.5, <f>o = 1 rad. 
The observation length is T s = T a = 10 days with no gaps and a noiseless dataset 
(we have set Sh = 1). In the upper panel we show IT (including AM sidebands) in 
solid grey, and 2T neglecting the AM sidebands (as defined by equation [2H)l in solid 
black. The dotted box in the upper panel indicates the area magnified shown in the 
lower left panel where we can identify the amplitude modulation sidebands for this 
single frequency modulated sideband. The dotted box in the lower left panel indicates 
the magnified area shown in the lower right panel where we have included a third 
evaluation of 2T through direct computation of equation [12] with frequency resolution 
1/T S . 



which are true for T s ^> P and at frequencies at or close to /„. We can now construct 
the ^-"-statistic via substitution of equations [151 and [161 into equation IT21 to obtain 

rp m 

2T ~ 2#" + A " ] + + A ' ] + 2C ^ A " + AsA ^ E J 2 n(Z)\W(f)\ 2 . (20) 

^ n=—m 

We have also chosen to neglect the amplitude modulation (AM) sidebands associated 
with each FM sideband for reasons of simplicity. As can be seen from the example 
jF-statistic spectrum shown in figured], there is significant additional jF-statistic located 
at the 8 AM sidebands per FM sideband, each separated by 1/ (sidereal day) Hz, and 
practical applications of this search technique should take advantage of this. 

4. Searching for sources with unknown frequency - the sideband search 

The search for continuously emitting sources is generally made more difficult when the 
intrinsic spin frequency of the source is unknown, and this is the case for LMXBs. We 
shall now describe a method for searching for these sources based on the incoherent 
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summation of sideband ^-"-statistic. This technique is a fast way to search broad 
frequency bands for sources with known sky position and orbital period. 



4-1. The sideband detection statistic 

The majority of jF-statistic values from a source in a binary system are spread broadly 
across a relatively wide bandwidth (~ 4irf aT s /P frequency bins) but are actually only 
present significantly in ~ 47r/oa isolated regions, each separated by 1/P Hz. We will 
define 

m 

<?(/)= £ *(/-/») ( 21 ) 

n=—m 

as an approximate search template where f' n is the frequency of the closest frequency 
bin to the frequency of the nth sideband, equal to f n = fo + n/P. We have therefore 
constructed a comb of unit amplitude "spikes" each separated by ~ 1/P Hz. We then 
convolve this template with the output of a wideband (sky position demodulated) T- 
statistic search, 2J r (/jj| giving us 

C(/) = 2jF(/)<g>g(/). (22) 

We should note that the number of "spikes" within the template is a function of the 
orbital radius and that an optimal value of C(f) will be obtained if using its true value. 
However, the power of this detection statistic is sensitive only to large ~ 50% mismatches 
between template orbital radius and the true value and is therefore not required as an 
exactly known parameter. This procedure is many orders of magnitude faster than the 
equivalent coherent matched filtered search but, as we shall describe in the following 
section, this increase in speed comes at the cost of reduced sensitivity. An example of 
the sideband statistic C(f) is shown in figure El 



4-2. The sideband search sensitivity 

The detection statistic C(f) (equation 122]) is simply the incoherent summation of 
sideband ^-"-statistics. In an idealised case of data free from gaps and with constant noise 
amplitude, our window function W(f) can be expressed analytically and the detection 
statistic becomes 

2 n t^ m I br(A/ n )] J 

We have used Af n = f n — f' n to represent the sideband frequency mismatches due to our 
finite frequency resolution. In practise Af Q will be drawn from a uniform distribution 
between the limits [—5f/2 — 5f/2), and for the purposes of this analysis we shall assume 
that this is the case for Af n . For the noise only case, 2T n is a ^-distributed random 
variable with mean and variance equal to 4 and 8 respectively. The simple summation 

|| In practise this is most efficiently evaluated via the convolution theorem. 
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of M = 2m + 1 such independent random variables results in a xlju'distributed random 
variable. The mean and variance of C for noise alone are then 

Mo(C) = 4M, (24) 
a 2 (C) = 8M, (25) 

respectively. When a signal is present in the noise, 2T n becomes a non-centrally \\~ 
distributed with a non-centrality parameter A n as a function of the random variable 
Af n . In order to approximate the distribution of C in this case we use 

2T n [A n (A/ n )] « |^(A/„)| 2 2^„(A n ), (26) 

where we have assumed that for values of |iy(A/ n )| 2 close to unity the change in the 
product of the T statistic with the window function is proportional to the relative 
change in the distribution mean. From this approximation it follows that the mean and 
variance for signal plus noise are 

pn(C) ^ p 2 opt (\W\ 2 ) + AM, (27) 

•ao. (# + ^ + -)( ^;aS; 16M2 ' 

-8p 2 opt (\W\ 2 )-16M, (28) 
respectively, where 

E A. = p 2 opt = F h 2 (t)W(t)dt. (29) 

and (...) indicates an averaging over the range of possible values of A/ n . By the 
central limit theorem we use the mean and variance approximations, in the case of 
signal plus noise, to describe the distributions of C as Gaussian. Figure. [2] shows the 
relationship between search sensitivity, parameterised as h \/T/ Sh, as a function of m 
and for various choices of frequency resolution 5f. We can also use Equations I2~4"l [23| I2T1 
and [28] to provide a rough guide to the sensitivity by approximating the SNR achievable 
with this search as 



^(C)- M o(C) . j / M 



For M ^> p 2 pt we see for a given sideband search SNR that 



ho oc M^y/Sh/T. (31) 

We should note that although this is an incoherent search technique, we still attain an 
ho sensitivity scaling proportional to T _1//2 , unlike most other incoherent searches. This 
is due to the fact that M ~ 2(27r/ofl) + 1 is a constant for a given source and therefore 
the M 1 / 4 factor in equation [31] remains constant for any choice of observation time. The 
clear comparison here is to the stack-slide search [6] where the statistical behaviour of 
the detection statistic is very similar to the behaviour of the sideband statistic. For 
stack-slide however, the analogue of the parameter M is the number of stacks, which 
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Figure 2. The sideband search sensitivity parameterised by h^^jT/ Sh for a 1% false 
alarm rate and 10% false dismissal rate averaged over the sky and as a function of 
m. We show in solid black, dark grey, and light grey, the curves corresponding to 
frequency resolutions of 1/T, 1/2T and 1/4T respectively. The dashed black curve 
represents our approximation to the solid black curve using Gaussian distributions to 
model noise and signal plus noise distributions. Note that for large m we see that the 
parameterised sensitivity is oc m 1 / 4 . 

have typically predefined length making the number of stacks proportional to the total 
observation time and consequently a T _1//4 scaling is achieved. 

5. A Follow up exploration of the jF-statistic likelihood 

We propose using an MCMC approach with our jF-statistic signal model components 
(equations [15] and [T6l) to explore the associated likelihood space. This is a natural 
extension to the sideband search in which the jF-statistic and its components have 
already been computed from our GW dataset. By selecting values of the previous 
"sideband" search output, C(f), above a predefined threshold, we identify isolated 
subspaces in our broadest parameter space dimension. These subspaces correspond 
to signal candidate frequencies and it is over this subset of frequencies that we perform 
the MCMC. 

5.1. An MCMC exploration of the J- -statistic likelihood 

In order to perform an MCMC we must calculate the likelihood. Let us first re- 
parameterise each nth complex ^-statistic components Fa and F^ n ' as 

F {n) = F (n) + - F (n) ^ p in) = F (n) + ^(n) ? (32) 

where n indexes the FM sidebands (— m < n < m). For data containing signals in the 
presence of Gaussian noise the real variables are Gaussian distributed with means 
given by 

fS n) = J n {Z)T{AA x + CA 2 )/4, ff > = - J n (Z)T(AA 3 + CA 4 )/4, 
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ff } = J n (Z)T(CA 1 + BA 2 )/A, ff } = - J n (Z)T(CA 3 + BA 4 )/4, (33) 
and governed by a covariance matrix 

where indicates a null 2x2 matrix. Therefore, the log likelihood function of the 
vector of ^-"-statistic components becomes 

00 ^("^2 , ^t?(") tfWy*? 



£(F) oc - ^ B [(FS n) - FS n) ) 2 + (F< n) - F 

n=—m 

+ a[(F^ -F^ n) ) 2 + A(Fi n) -F 



-2C 



(«)\2 
4 J 

( F H _ f< n) )(F< n) - F^ n) ) + (F< n) - F< n) )(F£ n) - Ff) . (35) 



We use an MCMC to efficiently explore this likelihood surface and to generate marginal 
posterior probability density functions (PDFs) for each of our search parameters. By 
using values of C(f) above a given threshold we can think of this as an approximate 
identification of small volumes of parameter space Vi likely to be associated with large 
values of the likelihood (or log-likelihood). By now treating fo as a search parameter 
within the MCMC we use Vi as an approximation to our prior knowledge on that space 
by setting p(f £ Vi) = 0. All subspaces Vi are then considered to form a continuous 
space on frequency such that a positive frequency jump from within subspace Vi large 
enough to fall just outside its upper boundary will land in subspace Vi+i. This way a 
single MCMC can be run and yield easily interpretable PDFs on our search space as 
opposed to separate PDFs for each frequency candidate which would leave us unable to 
make global parameter space inferences. 

For our simple LMXB model we define the primary search parameters as A = 
{fo, ho, cos t, ip, 0o, a, t p }. The form of the signal as described in equations [T5] and [16] 
shows us that the majority of ^-"-statistic is localised in M sidebands at frequencies 
fo + n/P. In the simplest cases this allows us to compute the log likelihood using only 
4M data points (ie. each of the 4 jF-statistic components sampled at the closest discrete 
frequency bin f' n to each sideband). Data used between sidebands would add only an 
identical constant factor to each computation of the log-likelihood for all A thus giving 
no useful discriminating information. By not using inter-sideband data the burden 
of the computation of £ is therefore reduced by a factor of ~ T/P and is therefore 
independent of observation time. We should also note that for our simplistic phase 
model each computation of C requires only a recalculation of the jF-statistic component 
means (equation 1331) for each of the M sidebands. Figure. [3] shows a set of marginal 
posterior distributions for an example LMXB type problem. We should also note that 
this MCMC stage alone is suited to the search for GWs from AMXPs for which the 
frequency and orbital parameters are known apriori to relatively high accuracy. 
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Figure 3. The marginal posterior PDFs (solid black lines) for an example LMXB 
type problem and the corresponding sideband detection statistic (lower right). The 
broken grey lines in the PDF plots indicate the true parameter values. The detector 
noise is assumed known at = 1 and the observation time is 10 days which for the 
chosen signal parameters results in an optimal SNR of ?. The search band is 0.01 Hz 
from which the 20 largest values of C were used as candidates for the MCMC prior 
frequency ranges. 

6. Summary 

We have presented a new incoherent search technique designed for performing GW 
searches for known objects in binary systems with unknown or poorly known frequencies, 
specifically LMXBs. Accurate knowledge of the orbital period allows the incoherent 
summation of frequency modulated ^-"-statistic sidebands removing orbital phase as a 
search parameter and also relaxing the requirement for accurate knowledge of orbital 
radius. Of course, all of these benefits come at the price of reduced sensitivity to 
GW amplitude which we have shown to scale as ho oc (47r/ a) 1/4 T 1/2 implying that 
greatest sensitivity is achieved for short period (small orbital radius) binaries. Full 
implementation of the search pipeline will take advantage of the multi-interferometer 
JF-statistic [28J whereby prior to the summation of sideband power data sets from N 
independent detectors are combined coherently thereby increasing sensitivity by ~ \/N 
for little additional computational cost. We have also identified that the coherent follow 
up MCMC stage constitutes a complete search if applied to the AMXPs, a similar 
class of source to the LMXBs but for which the frequency and orbital parameters 
are well constrained apriori. We again stress that here we have applied the proposed 
pipeline only to a simple source model and neglected issues such as poorly known orbital 
period and sky position, non-zero eccentricity, realistic detector noise characteristics and 
more complex issues such as the detailed spin frequency evolution of accreting systems. 
However, we believe that solutions to these problems can be integrated into this search 
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pipeline and are described in detail in [THJ [TB|. Finally, we stress that although the 
final search stage is a coherent one, the ultimate sensitivity of the complete pipeline is 
limited by the threshold set in the incoherent "sideband" search stage. 
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